/*
 * datatransform.h
 *
 *  Created on: Mar 8, 2011
 *      Author: tqlong
 */

#ifndef DATATRANSFORM_H_
#define DATATRANSFORM_H_

template<typename DATA>
class Normalize
{
public:
  typedef DATA                    dataset_type;
  typedef typename DATA::vec_type vec_type;
public:
  vec_type mean_;
  vec_type std_;
  double max_std_;
  dataset_type* data_;
public:
  void operator()(dataset_type& data)
  {
    std::cerr << "Normalizing...\n";
    data_ = &data;
    computeMeanStd();
    for (int i = 0; i < data_->n(); i++) {
//      vec_type tmp = (data_->at(i) - mean_);   // x = (x-mean) / std
      data_->at(i) -= mean_;
    }
//    std::cerr << "mean = \n" << mean_ << "std = \n" << std_;
  }
  void computeMeanStd()
  {
    mean_ = vec_type(data_->dim());
    std_ = vec_type(data_->dim());
    mean_.zeros();
    std_.zeros();
    int n = data_->n();
    for (int i = 0; i < n; i++) {
      mean_ += data_->at(i)/n;
      std_ += (data_->at(i) % data_->at(i))/n;
    }
    std_ -= mean_ % mean_;
    max_std_ = 0.0;
    for (int i = 0; i < data_->dim(); i++) {
      std_[i] = sqrt(std_[i]);
      if (max_std_ < std_[i]) max_std_ = std_[i];
    }
    if (max_std_ < 1e-10) max_std_ = 1;
  }
};

template<typename DATA>
class PCA
{
public:
  typedef DATA                    dataset_type;
  typedef typename DATA::vec_type vec_type;
  typedef std::vector<vec_type>   vec_vec_type;
  typedef std::vector<double>     dvec_type;
public:
  dataset_type* data_;
  vec_vec_type eigenvector_;
  dvec_type eigenvalue_;
  int n_, dim_, h_;
public:
  void operator()(dataset_type& data)
  {
    std::cerr << "PCA...\n";
    Normalize<DATA> normalize;
    normalize(data);
    data_ = &data;
    n_ = data_->n();
    dim_ = data_->dim();
    h_ = data_->n();

    double total_energy = 0.0;
    for (int i = 0; i < n_; i++)
      for (int d = 0; d < dim_; d++)
        total_energy += (1.0/n_)*data_->at(i).get(d)*data_->at(i).get(d);

    eigenvector_ = vec_vec_type(n_);
    eigenvalue_ = dvec_type(n_);
    double current_energy = 0.0;
    for (int p = 0; p < n_; p++) {
      bool no_change = false;
      vec_type vp = vec_type(dim_);
      double lambda;
      vp.randn();
      while (!no_change) {
        vec_type old = vp;
        mulCovariance(vp);
        orthogonalize(vp, p);
        lambda = arma::norm(vp, 2);
        vp /= lambda;
        no_change = arma::norm(vp-old,2) < 1e-10;
      }
      eigenvector_[p] = vp;
      eigenvalue_[p] = lambda;
      current_energy += lambda;
      std::cerr << "eigen " << p << " lambda = " << lambda << " = "
          << current_energy / total_energy << "\n";// << vp;
      if (current_energy / total_energy >= 0.90) {
        h_ = p+1;
        break;
      }
    }

    for (int i = 0; i < n_; i++)
      data_->at(i) = newVec(data_->at(i));
  }
  void mulCovariance(vec_type& w)
  {
    vec_type v = w;
    w.zeros();
    for (int i = 0; i < n_; i++) {
      w += (1.0/n_*arma::dot(v, data_->at(i))) * data_->at(i);
    }
  }
  void orthogonalize(vec_type& w, int p)
  {
    for (int i = 0; i < p; i++) {
      double dprod = arma::dot(w, eigenvector_[i]);
      w -= dprod*eigenvector_[i];
    }
  }
  vec_type newVec(const vec_type& x)
  {
    vec_type xnew(h_);
    for (int i = 0; i < h_; i++)
      xnew.ref(i) = arma::dot(eigenvector_[i], x);
    return xnew;
  }
};

#endif /* DATATRANSFORM_H_ */
